home *** CD-ROM | disk | FTP | other *** search
- // This may look like C code, but it is really -*- C++ -*-
- /*
- ************************************************************************
- *
- * Linear Algebra Package
- *
- * The present package implements all the basic algorithms dealing
- * with vectors, matrices, matrix columns, etc.
- * Matrix is a basic object in the package; vectors, symmetric matrices,
- * etc. being considered as a matrix of a special type.
- *
- * Elements of the matrix are arranged in memory in the COLUMN-wise
- * fashion (in the FORTRAN spirit). In fact, it makes it very easy to
- * feed the matrices to FORTRAN procedures, which implement more
- * elaborate algorithm.
- *
- * Unless otherwise specified, matrix and vector indices always start
- * with 1, spanning up to the specified limit.
- *
- ************************************************************************
- */
-
- #pragma once
- #pragma interface
-
- #include "myenv.h"
- #include <std.h>
-
- #define REAL float // Scalar field of the Linear Vector
- // space
-
- typedef REAL (*USER_DEFINED_REAL_F)(REAL);
- typedef double (*USER_DEFINED_DOUBLE_F)(double);
-
- class Vector; // Vector over the real domain
-
- class MatrixRow; // A row of the matrix
- class MatrixColumn; // A column of the matrix
- class MatrixDiag; // The diagonal of the matrix
- class MatrixPivoting; // For the determinant evaluation
-
-
- /*
- *------------------------------------------------------------------------
- * Basic class of matrix
- */
-
- class Matrix // General type matrix
- {
- friend class MatrixRow;
- friend class MatrixColumn;
- friend class MatrixDiag;
- friend class MatrixPivoting;
-
- private: // Private part
- int valid_code; // Validation code
- const int MATRIX_val_code = 575767; // Matrix validation code
-
- protected: // May be used in derived classes
- short nrows; // No. of rows
- short ncols; // No. of columns
- short row_lwb; // Lower bound of the row index
- short col_lwb; // Lower bound of the col index
- char * name; // Name for the matrix
- int nelems; // Total no of elements, nrows*ncols
- REAL ** index; // index[i] = &values(0,i) (col index)
- REAL * elements; // Elements themselves
-
- void allocate(const short nrows, const short ncols,
- const short row_lwb = 1, const short col_lwb = 1);
-
- public: // Public interface
-
- // Constructors and destructors
- // Make a blank matrix
- Matrix(const short nrows, const short ncols); // Index start with 1
- Matrix(const short row_lwb, const short row_upb, // Or with low:upper
- const short col_lwb, const short col_upb); // boundary specif.
-
- Matrix(const Matrix& another); // Make a new blank matrix a la
- // the other one
- Matrix(const char * file_name); // Read the matrix from the file
-
- ~Matrix(); // Destructor
-
- void set_name(const char * name); // Set a new matrix name
-
- // Erase the old matrix and create a
- // new one according to new boundaries
- void resize_to(const short nrows, const short ncols); // Indexation @ 1
- void resize_to(const short row_lwb, const short row_upb, // Or with low:upper
- const short col_lwb, const short col_upb);// boundary specif.
- void resize_to(const Matrix& m); // Like other matrix m
-
-
- void is_valid() const
- { assure(valid_code == MATRIX_val_code,"Invalid matrix"); }
-
- // Status inquires
- int q_row_lwb() const { return row_lwb; }
- int q_row_upb() const { return nrows+row_lwb-1; }
- int q_nrows() const { return nrows; }
- int q_col_lwb() const { return col_lwb; }
- int q_col_upb() const { return ncols+col_lwb-1; }
- int q_ncols() const { return ncols; }
-
- int q_no_elems() const { return nelems; }
-
- const char * q_name() const { return name; }
-
- // Individual element manipulations
- REAL& operator () (const int rown, const int coln) const;
-
-
- // Element-wise matrix operations
-
- // Matrix-scalar arithmetics
- // Modify every element of the
- // Matrix according to the operation
- Matrix& operator = (const double val); // Assignment to all the elems
- Matrix& operator -= (const double val); // Add to elements
- Matrix& operator += (const double val); // Take from elements
- Matrix& operator *= (const double val); // Multiply elements by a val
-
- // Comparisons
- int operator == (const double val) const;
- int operator < (const double val) const;
- int operator > (const double val) const;
-
- // Other element-wise matrix operations
- Matrix& clear(void); // Clear the matrix (fill out with 0)
- Matrix& abs(void); // Take an absolute value of a matrix
- Matrix& sqr(void); // Square each element
- Matrix& sqrt(void); // Take a square root
-
- Matrix& user_function(USER_DEFINED_REAL_F uf); // Apply a user function
- Matrix& user_function(USER_DEFINED_DOUBLE_F uf); // to each element
-
- // Element-wise operations on two matrices
- Matrix& operator = (const Matrix& source); // Assignment
- // Arithmetics
- friend Matrix& operator += (Matrix& target, const Matrix& source);
- friend Matrix& operator -= (Matrix& target, const Matrix& source);
- friend Matrix& add(Matrix& target, const double scalar,const Matrix& source);
- friend Matrix& element_mult(Matrix& target, const Matrix& source);
- friend Matrix& element_div(Matrix& target, const Matrix& source);
-
- // Comparisons
- friend int operator == (const Matrix& im1, const Matrix& im2);
- friend void compare(const Matrix& im1, const Matrix& im2,
- const char * title);
- friend void are_compatible(const Matrix& im1, const Matrix& im2);
-
-
- // True matrix operations
- // (on a matrix as a whole)
-
- // Transposition
- Matrix transpose(void) const;
-
- Matrix& operator *= (const Matrix& source); // Inplace multiplication
- friend Matrix operator * (const Matrix& A, const Matrix& B);
-
- // Matrix& operator *= (const DiagMatrix& diag); // Multiply by diagonal matrix
-
-
-
- // Compute matrix norms
- double row_norm(void) const; // MAX{ SUM{ |M(i,j)|, over j}, over i}
- double norm_inf(void) const // Alias, shows the norm is induced
- { return row_norm(); } // by the vector inf-norm
- double col_norm(void) const; // MAX{ SUM{ |M(i,j)|, over i}, over j}
- double norm_1(void) const // Alias, shows the norm is induced
- { return col_norm(); } // by the vector 1-norm
- double e2_norm(void) const; // SUM{ m(i,j)^2 }, Note it's square
- // of the Frobenius rather than 2. norm
-
- friend double e2_norm(const Matrix& m1, const Matrix& m2);
- // e2_norm(m1-m2)
-
- double determinant(void) const; // Matrix must be square one
-
- // Make a matrices of special kind
- Matrix& unit_matrix(void); // Matrix needn't be a square
- Matrix& hilbert_matrix(void); // Hilb[i,j] = 1/(i+j-1); i,j=1..max
-
- // Row, Column operations
- MatrixRow a_row(const int rown) const; // Get a matrix row
- MatrixColumn a_column(const int coln) const; // Get a matrix column
- MatrixDiag diag(void) const; // Get a matrix diagonal
-
- // I/O: write, read, print
- // Write to a file
- // "| command name" is OK as a file
- // name
- void write(const char * file_name,const char * title = "") const;
- void info(void) const; // Print the info about the Matrix
- void print(const char * title) const; // Print the Matrix as a table
-
- };
-
- /*
- *------------------------------------------------------------------------
- * Inline Matrix operations
- */
-
- inline Matrix::Matrix(const short no_rows, const short no_cols)
- {
- allocate(no_rows,no_cols);
- }
-
- inline Matrix::Matrix(const short row_lwb, const short row_upb,
- const short col_lwb, const short col_upb)
- {
- allocate(row_upb-row_lwb+1,col_upb-col_lwb+1,row_lwb,col_lwb);
- }
-
- // Make a new Matrix like the
- inline Matrix::Matrix(const Matrix& old) // old one
- {
- old.is_valid();
- allocate(old.nrows,old.ncols,old.row_lwb,old.col_lwb);
- }
-
- // Resize the matrix to accomodate to a pattern
- inline void Matrix::resize_to(const Matrix& m)
- {
- resize_to(m.q_row_lwb(),m.q_row_upb(),m.q_col_lwb(),m.q_col_upb());
- }
-
- inline REAL& Matrix::operator () (const int rown, const int coln) const
- {
- is_valid();
- register int arown = rown-row_lwb; // Effective indices
- register int acoln = coln-col_lwb;
-
- if( arown >= nrows || arown < 0 )
- _error("Row index %d is out of Matrix boundaries [%d,%d]",
- rown,row_lwb,nrows+row_lwb-1);
- if( acoln >= ncols || acoln < 0 )
- _error("Col index %d is out of Matrix boundaries [%d,%d]",
- coln,col_lwb,ncols+col_lwb-1);
-
- return (index[acoln])[arown];
- }
-
- inline Matrix& Matrix::clear(void) // Clean the Matrix
- {
- is_valid();
- memset(elements,0,nelems*sizeof(REAL));
- return *this;
- }
-
- inline void are_compatible(const Matrix& im1, const Matrix& im2)
- {
- im1.is_valid();
- im2.is_valid();
-
- if( im1.nrows != im2.nrows || im1.ncols != im2.ncols ||
- im1.row_lwb != im2.row_lwb || im1.col_lwb != im2.col_lwb )
- im1.info(), im2.info(), _error("The matrices above are incompatible");
- }
-
- // Assignment
- inline Matrix& Matrix::operator = (const Matrix& source)
- {
- are_compatible(*this,source);
- memcpy(elements,source.elements,nelems*sizeof(REAL));
- return *this;
- }
-
- /*
- *------------------------------------------------------------------------
- * Friend classes - MatrixRow, MatrixCol, MatrixDiag
- */
-
- class MatrixColumn // Special representation of a Col of the
- { // matrix
- friend class Matrix;
- friend class Vector;
-
- const Matrix * matrix; // The matrix I am row of
- REAL * ptr; // Pointer to the a[0,i]
-
- // Private constructor
- // Note there are no public constructors;
- // MatrixColumns are always
- // created via the Matrix::a_col()
- // function
- MatrixColumn(const Matrix* mp, const int col);
-
- public:
- // Note how to construct MatrixColumn of
- // the matrix 'm':
- // m.a_column(10)
-
- ~MatrixColumn() {}
-
- // Assign a value to all the elements
- // of the Matrix Col
- friend void operator = (const MatrixColumn& mc, const int val);
- // Modify the elements in the col
- friend void operator += (const MatrixColumn& mc, const int val);
-
- void operator = (const Vector & vec); // Assign a vector to a matrix col
- };
-
- // Construct the MatrixColumn
- inline MatrixColumn::MatrixColumn(const Matrix * mp, const int col)
- {
- matrix = mp;
- matrix->is_valid();
-
- register int colind = col - matrix->col_lwb;
-
- if( colind >= matrix->ncols || colind < 0 )
- _error("Column #%d is not within the matrix",col,((*matrix).info(),0));
-
- MatrixColumn::ptr = &(matrix->index[colind][0]);
- }
-
- inline MatrixColumn Matrix::a_column(const int coln) const
- return c(this, coln)
- {}
-
-
- class MatrixRow // Special representation of a Row of the
- { // matrix
- friend class Matrix;
- friend class Vector;
-
- const Matrix * matrix; // The matrix I am row of
- const int inc; // if ptr=@a[row,i], then
- // ptr+inc = @a[row,i+1]
- // Since elements of a[] are stored
- // col after col, inc = nrows
- REAL * ptr; // Pointer to the a[row,0]
-
- // Private constructor
- // Note there are no public constructors;
- // MatrixRows are always
- // created via the Matrix::a_row()
- // function
- MatrixRow(const Matrix* mp, const int row);
-
- public:
- // Note how to construct MatrixRow of
- // the matrix 'm':
- // m.a_row(10)
-
- ~MatrixRow() {}
-
- // Assign a value to all the elements
- // of the Matrix Row
- friend void operator = (const MatrixRow& mr, const int val);
- // Modify the elements in the row
- friend void operator += (const MatrixRow& mr, const int val);
-
- void operator = (const Vector & vec); // Assign a vector to a matrix row
- };
-
- // Construct the MatrixRow
- inline MatrixRow::MatrixRow(const Matrix * mp, const int row)
- : inc(mp->nrows)
- {
- matrix = mp;
- matrix->is_valid();
-
- register int rowind = row - matrix->row_lwb;
-
- if( rowind >= matrix->nrows || rowind < 0 )
- _error("Row #%d is not within the matrix",row,((*matrix).info(),0));
-
- MatrixRow::ptr = &(matrix->index[0][rowind]);
- }
-
- inline MatrixRow Matrix::a_row(const int rown) const
- return c(this, rown)
- {}
-
-
- /*
- *------------------------------------------------------------------------
- * Vector as a n*1 matrix
- */
-
- class Vector
- : public Matrix
- {
- friend class MatrixColumn;
- friend class MatrixDiag;
-
- public:
- Vector(const short n); // Specify a blank vector for a given
- // dimension, with indexation at 1
- Vector(const short lwb, const short upb); // Specify a general lwb:upb vector
- Vector(const Vector& another); // Make a vector by example
-
- // Make a vector and assign init vals
- Vector(const short lwb, const short upb, // lower and upper bounds
- double iv1, ... // DOUBLE numbers. The last arg of
- ); // the list must be string "END"
- // E.g: Vector f(1,2,0.0,1.5,"END");
-
- ~Vector() {}
-
- // Erase the old vector and create a
- // new one according to new boundaries
- void resize_to(const short n); // Indexation @ 1
- void resize_to(const short lwb, const short upb); // lwb:upb specif
- void resize_to(const Vector& v); // like other vector
-
- REAL & operator () (const int index) const;
-
- // Listed below are specific vector operations
- // (unlike n*1 matrices)
-
- // Status inquires
- int q_lwb(void) const { return row_lwb; }
- int q_upb(void) const { return nrows + row_lwb - 1; }
-
- // Compute the scalar product
- friend double operator * (const Vector& v1, const Vector& v2);
-
- // Vector norms
- double norm_1(void) const; // SUM{ |v[i]| }
- double norm_2_sqr(void) const; // SUM{ v[i]^2 }
- double norm_inf(void) const; // MAX{ |v[i]| }
-
- Vector& operator = (const double val)
- { Matrix::operator =(val); return *this; }
- };
-
- // Create a blank vector of a given size
- inline Vector::Vector(const short n)
- : Matrix(n,1)
- {}
-
- // Create a general blank vector
- inline Vector::Vector(const short lwb, const short upb)
- : Matrix(lwb,upb,1,1)
- {}
-
- // Make a vector by example
- inline Vector::Vector(const Vector& another)
- : Matrix(another.q_lwb(),another.q_upb(),1,1)
- {}
-
- // Erase the old vector body and create a
- // new one to accomodate a specified no.
- // of elements
- inline void Vector::resize_to(const short n)
- {
- Matrix::resize_to(n,1);
- }
-
- inline void Vector::resize_to(const short lwb, const short upb)
- {
- Matrix::resize_to(lwb,upb,1,1);
- }
-
- inline void Vector::resize_to(const Vector& v)
- {
- Matrix::resize_to(v.q_lwb(),v.q_upb(),1,1);
- }
-
- // Get access to a vector element
- inline REAL& Vector::operator () (const int ind) const
- {
- is_valid();
- register int aind = ind - row_lwb;
- if( aind >= nelems || aind < 0 )
- _error("Requested element %d is out of Vector boundaries [%d,%d]",
- ind,row_lwb,nrows+row_lwb-1);
-
- return elements[aind];
- }
-